Astronomy & Astrophysics manuscript no. pert Bondi-AA 


February 5, 2008 


(DOI: will be inserted by hand later) 





Perturbations of self-similar Bondi accretion 

Jose Gaite 

Institute de Matematicas y Fisica Fundamental, CSIC, Serrano 1 13bis, 28006 Madrid, Spain 
Received: 24 August 2005 / Accepted: 12 November 2005 

Abstract. The question of stability of steady spherical accretion has been studied for many years and, recently, the concept 
of spatial instability has been introduced. Here we study perturbations of steady spherical accretion flows (Bondi solutions), 
restricting ourselves to the case of self-similar flow, as a case amenable to analytic treatment and with physical interest. We fur- 
ther restrict ourselves to its acoustic perturbations. The radial perturbation equation can be solved in terms of Bessel functions. 
We study the formulation of adequate boundary conditions and decide for no matter-flux-perturbation conditions (at the Bondi 
radius and at r = 0). We also consider the problem of initial conditions and time evolution, in particular, of radial perturbations. 
No spatial instability at r = is found. The time evolution is such that perturbations eventually become ergodic-like and show 
no trace of instability or of acquiring any remarkable pattern. 
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1. Introduction 

The Bondi solutions for stationary spherical accretion onto a 
compact object (Bondi, 1952) have been the basis for many 
, analytical (or semi-analytical) studies of spherical accretion. 

Since there is a set of solutions, it has been a moot point how 
. to select the most adequate among them. Bondi suggested that 
' stability criteria would be useful for this purpose. However, 
neither analytical studies (Garlick, 1979; Petterson, SiUc & 
Ostriker, 1980) nor numerical studies (RufFert, 1994) have de- 
, tected any instatibility, in the sense that linear modes have am- 
■ plitudes increasing with time. Garlick's argument is very gen- 
, eral and powerful because it is based on energy conservation. 
' Kovalenko & Eremin (1998) pointed out another possible type 
of instability, namely, in space rather than in time. They stated 
the nature of this instability as given by a perturbation whose 
amplitude relative to the unperturbed quantity grows infinitely 
with decreasing radius, so that one expects that it reaches the 
nonlinear regime. Apart from this spatial notion of stability, 
Garlick (1979) had proposed before "to take a more sophisti- 
cated approach to the stability problem" (in the time domain): 
the perturbation, independently of its global behaviour, could 
be increasing in a small region (anywhere along the flow). 
However, he ruled this out as "unphysical". Garlick (1979) also 
considered briefly the initial value problem. 

Our purpose in this paper is to pay special attention to the 
problem of time evolution of perturbations of Bondi accretion 
and to go as far as we can with analytic methods. Given that the 
general problem of initial conditions and time evolution is in- 
trinsically difficult, we are compelled to make a number of sim- 
plifications. First of all, we will consider a particularly simple 
form of Bondi accretion; self-similar Bondi accretion. It occurs 



for the value of the poly tropic index y = 5/3, a limit case with 
physical interest (or it holds approximately as an intermediate 
asymptotic limit for y < 5/3). As an additional simplification, 
we will restrict ourselves to acoustic perturbations, neglecting 
entropic perturbations and vorticity. This is consistent, as gen- 
eration of entropy and vorticity are related, and they are simply 
advected with the flow and eventually vanish (Garlick, 1979). 
Moreover, to carry out a complete study of time evolution we 
will restrict ourselves to radial perturbations. In this case, the 
analytic treatment is simpler, as is the graphical representation 
of the results. 

We remark that a possible role of entropy and vortic- 
ity perturbations has been explored recently (Foglizzo, 2001). 
However, in accord with the traditional view, no instabilities 
can arise in polytropic flow in absence of an external source of 
entropy or vorticity. We consider no external sources. 

We divide the paper in three main sections: The first one 
(numbered second section) is a review of Bondi accretion and, 
furthermore, it introduces the self-similar solutions. The next 
section is devoted to the study of perturbations, in particu- 
lar, acoustic modes, and to formulate the appropriate bound- 
ary problem of differential equations. We follow previous work 
but we try to be more systematic. Unfortunately, the bound- 
ary problem does not lead to a Sturm-Liouville problem, even 
though is connected with one. The last main section focuses 
on radial perturbations, like the papers of Petterson, Silk & 
Ostriker (1980) and of Theuns & David (1992). In fliat section, 
we precisely formulate and solve the corresponding eigenvalue 
problem of ordinary differential equations (ODE's), finding an 
orthogonal system of eigenfunctions. Hence we solve the prob- 
lem of initial conditions. We end by discussing the relevance of 
our results for the problem of stabifity of Bondi accretion. An 
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appendix exposes the asymptotic form of small wave-length 
perturbations (WKB approximation). 

2. Steady spherical accretion: the Bondi solutions 

The problem of steady spherical accretion is suitable for an- 
alytical treatment and is the basis for more complicated and 
reaUstic problems (Frank, King & Raine, 1992; Shore, 1992). 
A compact object (e.g., a star) of mass M is at rest at the ori- 
gin (r = 0). An infinite cloud of gas, which is at rest at infinity 
and with density poo and pressure poo, is accreted by the com- 
pact object. The motion of the gas is assumed to be spherically 
symmetric and steady, so the fluid partial difl'erential equa- 
tions (PDE's) become ODE's in the radius r. Furthermore, the 
change of mass of the compact object, the efi'ects of radiation 
and viscosity are neglected. Then the equations describing the 
gas are 

d T 

— (rV) = , 



dv 



dr 

I dp GM 



(1) 

(2) 



dr p dr r^ 
with p the gas density field, v the gas velocity field, p the gas 
pressure and G Newton's gravitational constant. The continu- 
ity equation (1) can be integrated immediately. We must add 
an equation of state of the fluid, which we choose polytropic, 
following Bondi (1952). Then we have a set of three equations: 

4nr^pv = %i , (3) 

dv GM 1 dp 

V— = 2 -J- ' (4) 

dr r^ p dr 

P = 'Kip'' . (5) 

with TCi and "7^2 constants (characterizing the flow), and the 
polytropic exponent 1 < y < 5/3. With a suitable choice of 
7, the polytropic equation is equivalent to assuming that there 
is no heat transfer (adiabatic flow). The speed of sound, to be 
used later, is given by = dp I dp - 'Kijp'^'^. 

Integrating the Euler equation (4) over r one gets 
Bemouilli's equation. Furthermore, the polytropic equation (5) 
allows one to integrate dp, so the problem boils down to two al- 
gebraic equations, which can be expressed in non-dimensional 
form by using the sound velocity in the gas at infinity Coo. Let 

us define the non-dimensional variables 
r 

(6) 

P 



y = 



poo 
V 



u 



A = 



(7) 
(8) 

; (9) 

A7fP?-pooCoo 

that is, r is measured in units of the accretion (or Bondi) radius, 
H - GMjcl^, V in units of the speed of sound Coo, and p in units 
of Poo- The non-dimensional form of the equations of motion 
is: 

x^yu = X, (10) 
yT-^ - 1 1 



m2 

h 

2 y 



Solving for u in the first equation and substituting in the second, 
we obtain 



yx-v'+«(y'"-i) = x 



-1 



(12) 



where we have introduced the adiabatic index n - 1 /(y- 1). As 
an equation for x, it is a quartic algebraic equation, so it admits 
an algebraic expression for its solutions (however complex, see 
Theuns & David, 1992). 

2.1. Self-similar Bondi solutions 

The preceding algebraic relation between y and x (12) can be 
simplified in some cases. Let us assume first that the radial dis- 
tance is small, namely, that x"' » 1, so that the term n in 
Eq. (12) can be neglected (in comparison with x '). The condi- 
tion that the density is large, i.e., >■ » 1 , is equivalent (assuming 
that y is not too close to one). Then we have 

A2 



—x-^y-^+ny^'"x 
2 



1 



(13) 



For convenience, we define the variable z = y^^"x and rewrite 
this equation as 

yA;-3+2V^"+nz= 1 . (14) 

This is still a complicated relation between x and z, so we must 
look for further simplification. 

We see that in the particular case n = 3/2, Eq. (14) is 
just an algebraic equation in the sole variable z, with a solu- 
tion that only depends on the non-dimensional accretion rate 
A. Therefore, in this case, the non-dimensional density takes 
the form y - z(A)^^^ x^^^^ and the non-dimensional velocity 
takes the form u - AziAy^^^ x^^^^. We see that both density 
and velocity are power laws of the radius. More complicated 
solutions of Eq. (14) can be considered but we will focus on 
the case n = 3/2 ^ y = 5/3, which provides simple solutions 
and is realistic since it corresponds to adiabatic flow of perfect 
monoatomic gases. 

Then we write 

-3/2 



p(r) - ar 
v{r) = ^r-^l^ , 
with a and jS such that 

In this case, the speed of sound is 

c\r) = ^2 ^ = cr^r-' , 

so the Mach number is given by 



(15) 
(16) 

(17) 
(18) 
(19) 



1 



(11) 



c(r) cr 

which is constant. So each solution is charaterized by either 
{a,fi,a}or{'Ku'K2,M}. 

We remark that self-similar Bondi solutions have no dis- 
continuities and are either subsonic or supersonic everywhere. 
On account of the boundary condition at infinity, it is natural to 
discard the supersonic self-similar solutions. 
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2.1 .1 . Intermediate asymptotic limit when y < 5/3 

It is commonly observed that nonlinear equations that do not 
have similarity solutions at the outset develop them in an inter- 
mediate asymptotic regime, that is, a regime between two very 
different scales where both scales can be neglected (Barenblatt, 
1996). In the problem of steady spherical accretion by a com- 
pact object, we have two fundamental scales, namely, the Bondi 
radius K and the sonic radius. The latter is defined by the im- 
plicit equation 



GM 
2c2(r,) 



(20) 



(Frank, King & Raine, 1992; Shore, 1992). This equation can 
be written in non-dimensional variables as 



Zs(.Xs) — ys(.Xs) ^ Xs — , 



1 



(21) 



equation of state. In particular, they may be adiabatic perturba- 
tions; we shall henceforth denote thermodynamic functions as 
if this is the case. Then the equations for Unear perturbations 

are 



-iSpv+pdv) = , 

at 

^ +iv-V)6v + i6v-V)v = -V6x , 

ot 



where 



(23) 
(24) 

(25) 



is the enthalpy perturbation (the speed of sound is given by 
= TCzyp^"^). Therefore, we have a system of four linear 
PDE's. 



where >'.v(-«.s) is the solution of Eq. (12) corresponding to 
Writing Eq. (12) in terms of Xj and Zj = 1 /2, we have 



22«-1^2 ^-3+2n 



nxs + I - - 



(22) 



If2>n>3/2<=>3/2<7<5/3 this equation has zero, one 
or two solutions, dependig on the value of the non-dimensional 
accretion rate A (given n). As A increases, the number of so- 
lutions increases as well, being one at the critical value Amax, 
which is only function of n. This critical value is called Amax 
because it is the maximal value such that the flow is continuous 
(Theuns & David, 1992). The value 

_ 1 3 
-*^.s('^max) - ^ ~ 

2 An 

is only function of n and approaches zero as n ^ 3/2. 

If n > 3/2 but not much larger (y < 5/3 but not much 
smaller), and A is not much smaller than Amax, the solution of 
Eq. (14) for z will have a dependence on x but only logarithmic, 
such that the power laws (15) and (16) will hold with correc- 
tions proportional to log(r/r5). In other words, in the intermedi- 
ate asymptotics « r « ??, the power laws are approximately 
correct (so the Mach number M hardly grows with the flow). 

For definiteness, we will set y = 5/3 in the following, keep- 
ing in mind that the results are of somewhat wider applicability. 

3. Linear perturbations of the self-similar Bondi 
solutions 

Let us assume that the steady sphericaUy symmetric solutions 
are perturbed, such that the total velocity and density fields be- 
come p(r)-i-dp(jir, t) saAv{r)+6v{x, t), with p(r) and v(r) = v{r)f 
the self-similar solutions, given by Eqs. (15) and (16) in the 
self- similar case (r is the unit radial vector). 

The independent perturbations are 6p, Sv and the "entropy 
perturbation" 6s, where the "entropy" is defined as 



log- 



yp 

(Kovalenko and Eremin, 1998). We shall assume that the per- 
turbations are such that ds = 0, so that they keep the polytropic 



3.1. Solution for acoustic modes 

For the acoustic modes, 6s = and the perturbation in the 
velocity field can be written as 6v = V0. In this case, and given 
that Bondi flows are curl free, we can write Eq. (24) as the 
following differential equation for the scalar potential field (f> 



d(h 

+ V(v • V^) = -V6x . 

at 



(26) 



This equation can be integrated to obtain (Kovalenko & 
Eremin, 1998) 



Tt 



-i-v • V0 



(27) 



which provides 6x and hence 6p in terms of 0. Equation (23) 
can then be written in terms of as follows 



9 A 
-+vV + 
at p 



-V • (pV0) 



(28) 



[Eq. (24) of Kovalenko & Eremin (1998)]. This is an equation 
for non-homogeneous wave propagation. 

One tries a solution of equation (28) for the scalar potential 
^ by separation of the variables in spherical coordinates, that 
is, a solution in the (complex) form 



<l>(r,e,<f,,t)=R(f)Yije,ip)e- 



(29) 



Given this form of the scalar potential, the velocity perturbation 
is given by 

6V = {SVr, SVg, SV^) , 

with 



6vr(r, 0, <p, t) 
Svgir, 6, (f, t) 



Re[e 



dr 



Re[e-"''—dgYUe,ip)], 



6v^{r,e,(p,t) = Re[e 



t R(r) 



(30) 
(31) 
(32) 
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(after taking the real part of the complex quantities). Note that 
the angular components form a poloidal vector field and the 
form of the radial component ensures the absence of vorticity 
(Garlick, 1979). 

The differential equation satisfied by /?(r) is 

(cr2 -p^y^R" + ( + 2i/3a>r^'^)rR' + 



[-/(/ + l)cr^ + a/P + i/3cor^'^] R = 0. 



(33) 



Following standard methods of the theory of ODE's, the second 
order equation 



R"(r) + P{r)R'{r) + Q{r)R{r) = , 
with 

^^'^ = • 

can be simplified by the following change of variable: 

l„g«W = l„g«,)-iJ*PW, 

which integrated yields 



- ..-1/4 -2iPojr 



3/2 



R{r) = r"" exp 



3(cr2-,82) 



h(r). 



(34) 

(35) 
(36) 

(37) 

(38) 



In terms of this new dependent variable hir), the differential 
equation (33) becomes 



h"ir) + (air-^ + airMr) = , 
with 

3 /(/+I)cr2 



ai 



ai = 



16 o-^-p^ 

\2 



,0-2-/32 



(39) 



(40) 



(41) 



The preceding differential equation (39) has one regular 
singular point at r = and one irregular singular point at r = cx). 
In fact, it can be transformed into the Bessel equation by mak- 
ing r^''2 the dependent variable. So we obtain 

R(r) = r"^e-"''"\CiMKP'^) + C2NMP'^)] , (42) 



where 



2crcjj 



_ V[l + 16/(/+ 1)] 0-2-/32 

6Vo-2-y82 

and Ci , C2 are, in general, two arbitrary complex numbers. 

Note that the order v of the Bessel function in solution (42) 
does not depend on j3 and cr separately but on its ratio, that is, 
the Mach number. At = j6/o- (in addition to depending on the 



angular number /). If the values of M and / are generic such 
that V is non-integer, we can write the solution of the equation 
as 

R(r) = r^''^e-""'"\CiJyiKp'^) + C2J-y(Kr^'^)] ■ (43) 

In this form, the most general (complex) solution of the scalar 
potential is 



{CiJyiKr^'^) + C2J-yiKr^'^)). 



(44) 



Since we assume At < 1, the signs of fi and k are positive and 
V is real. 

Boundary conditions are needed in order to determine both 
Ci and C2 and the a>-spectrum. Notice that the constants C\ 
and C2 need to be specified for each solution, that is, they are 
indexed by /, m and the w-spectrum index. 

3.2. Boundary conditions 

The acoustic perturbations satisfy the wave equation (28), so a 
suitable boundary condition on a surface is needed. After the 
separation of variables in spherical coordinates, one needs to 
impose two boundary conditions on the radial differential equa- 
tion, as this equation is of second order (Morse & Feshbach, 
1953). So one assumes that the boundary conditions adapt to 
the spherical coordinates and, in particular, the boundaries are 
two spherical surfaces of inner radius ri and outer radius r2, re- 
spectively. The natural inner radius is either the accretor or the 
sonic radius. The natural candidate for outer radius is the Bondi 
(or accretion) radius K: as r gets larger than the thermal en- 
ergy of the nearly homogeneous gas becomes dominant over 
its gravitational energy and we assume that this homogeneous 
gas remains unperturbed. 

The choice of boundary conditions of Petterson, Silk & 
Ostriker (1980) (for radial perturbations) is that the flux per- 
turbation Sf should vanish at both ri and since the flux of 
matter through a sphere of radius r is constant in the Bondi so- 
lutions [Eq. (3)], it seems natural to impose that the perturba- 
tions do not change it on the inner and outer boundary. Unlike 
for radial modes, for general acoustic modes the variation of 
the matter flux current j = pv, namely, Sj = p6v + 6pv , de- 
pends on the angles 6, (f. The perturbed flux of matter through 
a sphere of radius r is given by the angular integral of its radial 
component, proportional to J dClYim{6,(p). This angular inte- 
gral vanishes if / ^ 0, so this condition is only useful for radial 
modes. 

Instead, Kovalenko & Eremin (1998) derive their boundary 
conditions from the law of conservation of acoustic energy 

^=V-W, (45) 
where the acoustic energy density is 

\2 



E = 



P 

2c2 



and the acoustic energy flux current is 



c2 dt 



V ^ - c^V0 -I- V • (v • V0) 

ot 



-dt(f>6j. 



(46) 



(47) 
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The acoustic energy flux through a spherical surface centred on 
the origin is given by 

J dQ.Wr = ^ J dQ. dt<p [v d,<p - c%<p + v^dr<p] . (48) 

This quantity is quadratic in (p and does not necessarily vanish 
for any / upon performing the angular integral. 

We see from Eq. (48) that sufficient conditions for the van- 
ishing of the acoustic energy flux are either 



dt(f> = 



or 



r^r-6j = — J- [v 5,0 - c^drcfi + v^dr(f>] = . 



(49) 



(50) 



The latter condition implies, of course, the vanishing of its inte- 
gral 6f over any angular domain, so it is the appropriate gener- 
alization of the condition of Petterson, Silk & Ostriker (1980). 
Hence, we choose it as our boundary condition; namely, this 
quantity should vanish at both r\ and r2. It can be expressed 
only in terms of Rir) as 



2 



[icovR + (c^ - v^)R'] = . 



(51) 



The discussion above is general, that is, it is not restricted 
to self-similar solutions. For self-similar solutions, the accretor 
and sonic radii vanish, so we take r\ = 0; on the other hand, the 
Bondi radius 'R—^oo. However, we need a finite outer radius to 
have a discrete spectrum, so we take rz - H finite, keeping in 
mind that the actual spectrum becomes continuous in the hmit 

The radial equation (33) for R{f) is not self-adjoint [the sub- 
ject of adjointness of ODE's and the Sturm-Liouville bound- 
ary problem is treated, for example, by Morse & Feshbach 
(1953) and Coddington & Levinson (1955)]. However, the 
transformed equation (39) for h{r) is patently self-adjoint. The 
boundary condition (51) in terms of h{r) reads 



r-^l\-h + 4rh'ir)] = 



(52) 



(at r\ and rz). Therefore, the boundary problem for h(r) is of 
Sturm-Liouville type and the eigenvalues and eigenfunctions 
fulfill the corresponding properties; in particular, the eigen- 
functions are orthogonal (Morse & Feshbach, 1953). This is 
natural, because h(r) is a combination of Bessel functions (with 
a factor r'^^). In contrast, the boundary problem for R(r) is not 
of Sturm-Liouville type and the eigenfunctions need not be or- 
thogonal (note that the k dependence of /i in Eq. (42) spoils the 
orthogonahty inherent to the Bessel functions). Having an or- 
thogonal system of eigenfunctions is important for the problem 
of initial conditions, namely, for finding the coefficients of the 
expansion corresponding to the initial data. 

We now apply Eq. (52) at ri =0 and rz. For the former, we 
need to take the limit r — > of Eq. (52). Writing it in terms 
of Bessel functions, with the dependent variable z = kP^^, we 
have 

Umz-i/«(Ci[/v(z) + 6z/;(z)]+ 

z->0 



In this limit, it is sufficient to consider the first term of the 
Taylor expansion of the Bessel functions, namely. 



r(v + i) 



(1 + 6v)z 



y-U6 



(54) 



When V > 1/6 (/ > 0) we are led to imposing C2 = to remove 
negative exponents (divergent as r — > 0). If v - 1/6 (corre- 
sponding to the radial modes), we must take Ci = instead. 
The outer boundary condition at rz yields 



JyiZ2) + 6z2J'iZ2) = 0, 



(55) 



which provides the allowed eigenvalues k„, and hence the oj„. 

Finally, it is important to analyse the relative variations of 
the density and the velocity. Their behaviour in the limit r — > 
is: 



6Vr 



-1/4+3V/2 



V 

^ ^ ^ ~ ^1/2^ _ ^-l/4+3v/2 

V V r 

^ = \[icoR(r) - vR'ir)] ~ r"^^^''^ , 
p 



(56) 
(57) 
(58) 



C2[J-yiz) + 6z/ ,(z)]) = . 



(53) 



for V > 1/6; for v = 1/6, 6vr/v ~ 6p/p ~ r. All these relative 
variations vanish as r — > 0, so there are no spatial instabilities 
at the inner radius and Unear perturbation theory is sound. 



4. Radial perturbations 

The simplest acoustic modes are the radial perturbations, stud- 
ied by Petterson, Silk & Ostriker (1980) and Theuns & David 
(1992), for general Bondi flow. Indeed, a more thorough an- 
alytical treatment is also possible in the self-similar case if 
we restrict ourselves to radial perturbations. Petterson, Silk & 
Ostriker (1980) obtained general stabiUty results: the problem 
of standing waves with vanishing flux boundary conditions has 
only real eigenvalues of u), and radially travelUng waves can be 
shown not to grow (by an energy argument). Theuns & David 
(1992) derived an associated Sturm-Liouville boundary prob- 
lem by decomposing the flux perturbation ^/ as if it were com- 
posed of modulated waves with a position dependent phase ve- 
locity. This decomposition leads to a self-adjoint equation for 
the amphtude. Their procedure is analogous to the transforma- 
tion that led us from the radial Eq. (33) to the self-adjoint equa- 
tion (39). 

Although the associated Sturm-Liouville boundary prob- 
lem yields the eigenfunctions of the original boundary prob- 
lem, it does not provide us with an orthogonality relation for 
them, which may exist or not. Fortunately, a specific mathe- 
matical treatment of the first order equations for radial pertur- 
bations allows us to find a definite scalar product under which 
the eigenfunctions are orthogonal. This is exposed in appendix 
A. Next, we show the solution to the first order boundary prob- 
lem, with the orthogonality relation, and then we proceed to 
solve the problem of initial conditions to study the time evolu- 
tion. 
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4.1. Solution of the boundary problem 

The general linear perturbation equations (23) and (24) re- 
stricted to radial modes give 



0, 



0. 



(59) 



(60) 



-icodp + [?"^(5p V +p ^v)j 

-icj 6v + dr\ — 6p + v6v\ 
\P I 

We denote y = {6f,6g), where 6f = r^iSpv + pSv) is the 
perturbation of the flow (but for the factor 4jt) and 6g = 
{c^lp)6p + v6v - icocp [according to Eq. (60)]. Solution (43) 
with V = 1/6 and Ci = (because of the boundary condition at 
r = 0) gives 

6 fir) = exp {-iMKr^'^) r^'^ i Js/eixr^'^) , (61) 

6g{r) = exp [-iMKr^'^) ^ r^'"^ J-i/eiKr^'^) , (62) 

where we have used the identity J-i/eiz) + 6z7^j^g(z) = 
-^zJs/eiz) to simpUfy 6f, and we have divided both com- 
ponents by the constant total flux so that 5f becomes non- 
dimensional. Now we impose the boundary condition 6f - 
at r = 1 (we normaUze the radial variable by dividing it by 
r2 - "R). It determines that the posibles values of k are the ze- 
ros of the Bessel function J^jf, and so yields the eigenvalues k„. 
Redefining the time unit such that t at, we have the modes 



Sfnir, f) = exp 



Sgnir, t) = exp 



where 



^„^Mr3/2+^(l-M2);| 



a 



3/2x 



(63) 



(64) 



are the zeros of the Bessel function J^/e- 
Furthermore, one can reverse the sign of k„, obtaining the con- 
jugate eigenfunction );* = {6f*,6gl), up to a constant. Hence 
we denote = and y^n - y*„- An additional constant 
mode is necessary, namely, yo = (0, 1) (arising from the limit 
of Eqs. (61) and (62) as -> 0). It corresponds to a constant 
velocity potential. So the full set of modes is {.y„))^_oo- 

These modes are orthogonal with respect to the product of 
Eq. (A.10),thatis, 



(yn,ym) 



Jo 



1 



-^2 



cr 



p r'l\6f: 6g„ + 6gl 6U) + a ?'^5g, 5gJ^ dr cc 6„,„ . (65) 

(To have a non-dimensional scalar product, we have introduced 
the constant or.) We have checked numericaUy the orthogonal- 
ity relation (65) for particular values of At and for modes up to 
n ^ 100. 



4.2. Initial conditions 

Having solved the boundary problem in terms of orthogo- 
nal modes, we are in position to solve the problem of ini- 
tial conditions. As initial conditions, we must take two func- 
tions 6 fir, 0) and 6gir, 0) that satisfy the boundary conditions, 



namely, (5/(0,0) = 5/(1,0) = 0. Then we expand y = (6f,6g) 
in modes as 



yir) = ^ c„y„ir) 



(66) 



Since y_„ = y*, we require c_„ = c* to have real y. So we can 
also write 



y(.r) =co)'o + 2Re 



Lji=l 



Using the scalar product of Eq. (65), we can obtain the coeffi- 
cients in the expansion (66): 



Cn = 



{yn,y) 
{yn,yn) 



(67) 



Either (5/(0, r), 6gi0, r)) or the set of coefficients {c„)^q con- 
tain the information on both density and velocity perturbations 
and constitute alternate definitions of the initial conditions. 



4.3. Time evolution 

Given the set of coeflicients {c„)^q, the solution of the bound- 
ary problem in terms of the modes (63) and (64) can be written 



y{r,t) = cojo +2Re 



£c„e-''^"'j«('-) 



Ln=l 



(68) 



Therefore, y{t,r) adopts an expansion similar to the one of 
3^(0, r), Eq. (66), but with time-dependent coefficients 

^n(f) — 



Note that the norm 

CO 

{y,y} = \cof'{yo,yo} + 2^ |c„P<y„,y„). 



(69) 



n=l 



is invariant in time. We can express this norm in terms of 
the variables 6p and 6v or, alternatively, in terms of dt^ = 
-c~6p/p - vSv and dv: 



iy,y) 



pi 2 

= ^ [(5,0)2 + (c'-v')(<5v)2]fl!r 



(70) 



Comparison with the expression for the energy density (46) 
shows that the norm is just the energy (but for a constant fac- 
tor), so it is naturally invariant in time. Therefore, the mode 
decomposition of the norm (69) expresses the total energy as 
a weighted sum of the energies of all modes. The "energy 
spectrum" {|c„P}^q (invariant in time) characterises the initial 
conditions in combination with the initial phases of the coef- 
ficients. Suitable initial conditions must be normalizable, that 
is, must have finite energy. From the asymptotic form (B.IO), 
{yn,yn) 0^ 1/(1 + 6n), we deduce that we must demand that 
lim„^oo Icnf = (otherwise, the asymptotic form of the series 
(69) diverges like the harmonic series, at the least). 

If we define the initial conditions by the energy spectrum 
and the initial phases of the coeflicients, the forms of 6f and 
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6g crucially depend on the correlation between those phases. 
In general, the evolution given by Eq. (68) is such that the cor- 
relation between the phases decreases with time and is finally 
lost as / — > oo. To be precise, this happens unless the eigen- 
values co„ have a common ratio, in which case the solution is 
periodic. But this does not happen in our case. In other words, 
the type of motion that we find is such that it leads to an evo- 
lution that is ergodic in a restricted sense, that is, an evolution 
that conserves the energy spectrum but which otherwise ex- 
plores the available phase space. A typical state is given by 
taking the coefficient phases to be random. It could happen that 
either 6f(r, f) or 6g(r, f) became concentrated about a value of 
r at some t (and they could even diverge) but it is extremely 
unUkely. 

We have studied the evolution of particular radial perturba- 
tions of self-similar Bondi flow with initial conditions that are 
sufficiently smooth. They must have an energy spectrum that 
decays rapidly as n — > oo. As this condition is time-invariant, 
smoothness is preserved in time. Functions such that they are 
concentrated in a sort of bell shape have special interest. A 
Gaussian profile is not allowed because of the boundary con- 
ditions {6f vanishes at r = 0, 1); but we can take polynomial 
bell-shaped functions, for example, fo„(x) = [4 x{l - x)Y, n = 
1,2, .. . (which are succesively more concentrated). The time 
evolution of initial configuration 6f{r,Q) - b4(r),6g{r,0) - 
is plotted at several times in Fig. 1. This function gives rise 
to an energy spectrum that decays very rapidly, so the modes 
with n < 50 correspond to a relative error in its norm (total 
energy) smaller than 10"^. We observe that 6g initially evolves 
fast, reaching a non-negligible value at r = in a very short 
time. We also observe that 6f(r,0) is almost reproduced at 
f ^ 1.7 (but 6g{r,0) = does not seem to be reproduced). 
Nevertheless, there is no periodicity and the reason why 6f 
seems periodic is that the power spectrum is very concentrated 
in a small number of modes that are distributed with sufficient 
regularity. At any rate, the phase correlation among them is de- 
stroyed for longer times, as the following plots in Fig. 1 prove. 
The last couple of plots show a "typical" long-time state, with 
random phases. 

5. Discussion 

We have examined some of the most interesting questions re- 
garding the instability of Bondi accretion. The first question is 
about the spatial stability of perturbations, in the limit r — > 0. 
For our self-similar solutions, the meaning of this limit is dif- 
ferent from the meaning in Kovalenko & Eremin's (1998), be- 
cause similarity implies that the accretor and sonic radii vanish. 
The behaviour of perturbations in the limit r ^ depends on 
the type of inner boundary condition. We have shown that the 
natural boundary condition is the no-energy-flux perturbation 
condition, which is equivalent to the vanishing of either the ve- 
locity potential time derivative or the radial matter current vari- 
ation. We have focused on the latter condition as the natural 
generalization of the boundary conditions of Petterson, Silk & 
Ostriker (1980) for radial perturbations. Under that condition, 
the relative variations of the physical variables, velocity and 
density, are bounded in the limit r -> 0. Therefore, the spatial 




Fig. 1. Time evolution of initial configuration <5/(r, 0) = 
b4(r),Sg(r,0) = 0, with modes up to n = 50 (corresponding to 
a relative error in the norm < 10"''). We have taken M = 1/2. On the 
left-hand side are plots of 5f(r, t) and on the right-hand side are plots 
of Sg(r,f), for f = 0,0.1, 1.7, 100. The bottom plots correspond to a 
realization of the same energy spectrum with random phases. 

stability of perturbations is ensured. This conclusion disagrees 
with an assertion by Kovalenko & Eremin (1998), namely, that 
some perturbations are spatially unstable, for example, a non- 
radial subcritical mode. We attribute this difl'erence to their im- 
posing the no-energy-flux condition only under time average 
(which seems incomplete). 

The second question is about the time evolution of per- 
turbations. We have seen that the evolution is ergodic-like, 
namely, the perturbations make a random walk in the section 
of the total phase space with given energy spectrum, losing 
memory of the initial conditions (the initial phases of the coef- 
ficients c„). Regarding Garlick's proposal of perturbations that, 
independently of their global behaviour, could be increasing in 
a small region, we can now conclude that that can actually hap- 
pen, but it is extremely unlikely: indeed, we have shown that 
an initial condition concentrated in a small region is forgotten 
in the evolution and not encountered again in a random walk in 
the big phase space. Furthermore, if we consider the non-linear 
coupling among the modes, with the corresponding transfer of 
energy, we deduce that the energy spectrum is not really con- 
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served, so the only truly conserved quantity is the total energy. 
Therefore, in the long run, the energy spectrum thermaUzes and 
the evolution becomes fully ergodic. In addition, the presence 
of viscosity (or other dissipative processes) would also trans- 
form the energy of perturbations into thermal energy. We see 
that an energy concentration in a small region corresponds to 
a situation of entropy spontaneously decreasing; so we fully 
agree with Garlick that this possibility must be ruled out as un- 
physical. 

We have restricted ourselves to acoustic perturbations of 
self-similar Bondi accretion, but many results must be appUca- 
ble to non-self-similar Bondi accretion. For example, the stabil- 
ity and ergodic-like nature of the evolution of acoustic pertur- 
bations, ultimately due to the reality of the spectrum of eigen- 
values (i)„ and to its genericity (in the sense of not having a 
uniform distribution). 

Acknowledgements. I am grateful to Carmen Molina-Pari's for conver- 
sations. My work is supported by the "Ram6n y Cajal" program and 
by grant BFM2002-01014 of the Ministerio de Educacion y Ciencia. 

Appendix A: Eigenvalue probiem for radiai 
perturbations 

The perturbation equations (23) and (24) become for radial 
modes 

d,6p + [r^iSp V + p 6v)\ = , (A. 1) 

d,6v + dr^^6p + v6vj = 0. (A.2) 
After performing the Fourier transform for the time variable, 
-ioj 6p+^dr [r\6p v+p 6v)\ = , (A.3) 

(A.4) 

We define the vector x - (r^6p, 6v) and the matrix 



-io) Sv + dr\ — 6p + v6v\ = . 



A = 



v(r) r- p{r) 
v(r) 



such that detA - - Q. Then we can write the ODE's 
(59) and (60) as the vector equation 



d{A ■ x) 
dr 



ICOX. 



Its adjoint equation is 
dr 



(A.5) 



(A.6) 



where the adjoint matrix is, in this case, just the transpose 
of A, and we have written the unknown vector function as y for 
later convenience. The derivation of this adjoint equation holds 
as long as the boundary conditions are such that y* -Ax vanishes 
at the ends of the interval of integration that defines the scalar 
product of functions. 



According to the general theory of ODE's (Coddington & 
Levinson, 1955), the eigenfunctions of the original ODE and 
the eigenfunctions of its adjoint equation are orthogonal, that 



IS, 



y'n-XmdrocSn 



(A.7) 



and their respective eigenvalues are conjugate. But the eigen- 
values are in fact real in our case due to a property of the matrix 
A, namely, it is transformed into its adjoint by the interchange 
of the two components of vectors. To make this property for- 
mal, let us define the matrix 



1 

1 



U. Then = f/ • A • f/, so we 



(A.8) 



U = 

such that - U and 
obtain from Eq. (A.6) 

d(U ■%) ■ *Ti - 
A — — = ioj„U -yn. 

On the other hand, defining 



y„ = A ■ x„ = {r {6p V + p 6v), — 6p -n V(5v)„, 

P 

we can write Eq. (A.5) as 

. dyn . 
A - — = ia)„y„. 
dr 



(A.9) 



Therefore, we can identify v„ - U-yn, that is, we have an equiv- 
alence between the original equation and its adjoint. In conse- 
quence, we have that ll>„ = w*. This constitutes a proof of the 
reality of eigenvalues different from those given by Petterson, 
Silk & Ostriker (1980) or Theuns & David (1992). 

Furthermore, given the preceding identification, we can 
write the orthogonality relation (A.7) as the orthogonality of 
eigenfunctions x„ or, alternately, of eigenfunctions y„; we shall 
use the latter, namely. 



U-A- 



■yr 



, dr oc 6n 



(A. 10) 



Note that the matrix U ■ A ' is self-adjoint and positive and 
plays the role of a metric in the space of vectors y. 

As we have pointed out above, the adjoint equation is valid 
if y*-A X - >•*•>• = y' y^+>'^ y' = at the ends of the integration 
interval (using the components of y = (y^,y^)). Therefore, suf- 
ficient boundary conditions are given by either component y^ 
or vanishing at both ends. We choose the former to vanish, 
as explained in the main text. 

Appendix B: Asymptotic form of acoustic modes 

B.1. Connection with WKB solution 

When the wave-length is small compared with the typical di- 
mensions of the zone in which waves propagate, the WKB so- 
lution is a good approximation. The WKB solution of the dif- 
ferential equations for radial modes was obtained by Petterson, 
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Silk &Ostriker (1980): 
1 



Sfir) 



Spir) 



y/cv 



A2 exp J"-^ jj , 



V + c 
v-c 



■ exp 10) 



■ exp 



f dr 
10) I 

J v + c 



(B.l) 



(B.2) 



In this equation, the two functions that form the linear com- 
bination with (arbitrary) coefficients Ai and A2 represent the 
outgoing and incoming waves, respectively. Substituting the 
self-similar expressions of v and c and performing the integrals, 
we obtain the corresponding expressions of the WKB solutions 
Sfir) and 6pir). From them we can further obtain 



Sgir) = - r 
a 



-1/2 



Ai exp ; 



■ A2 exp ; 



(B.3) 



308 + £r) "3(j8-cr)J 

This formula is valid for large o) or, equivalently, large r. 

The two component waves are analogous to the functions 
Jy and A^v in the radial solution (42). To find the precise rela- 
tion between them, it is convenient to express the radial solu- 
tion (42) in terms of the complex combinations of the Bessel 
functions, namely, the Hankel functions H^y^ and H^y^ (Morse 
& Feshbach, 1953). The asymptotic forms of the Hankel func- 
tions for large values of their arguments is 



(B.4) 



(B.5) 



Expressing Eq. (42) in terms of Hankel functions (with two 
new constants C[ and C'^) and using the preceding asymptotic 
forms, we have 

Rir) ~ /''*e-'^"^"r-^'\C{ expimr^'^) + C'2expi-iKp'^)l (B.6) 



still with /i 



2/3m 



2aoj 



Note that this asymp- 



totic form becomes v-independent and therefore l-independent. 
Given that 



we have that 



K = - 



2/3aj 



2cra> 



2a> 



3(0-2 _ ^2) 3(^2 _ ^2) 3(J3±a-)' 



Rir) 



-1/2 



2ia)r^'^ ^, 2iwr3/2 
Cj exp — — -I- exp 



(B.7) 



3(/3 + (T) ' •^3(J3-(r)\ 

This last form is equivalent to the WKB form (B.3) but is valid 
for non-radial modes as well. 

B.2. Asymptotic form of radial eigenfunctions 

For n » 1 we can substitute in Eqs. (63) and (64) the asymp- 
totic forms of the Bessel functions for large values of their ar- 
guments (Morse & Feshbach, 1953): 

Jy{x) = ^cos - Y - ^) + Oix-^'^) 





Fig. B.l. Real parts of radial modes with n = I and n = 4 (solid lines) 
compared with their asymptotic trigonometric forms (dashed Unes). 
Note that the difference diminishes as r grows. 



[consistently with Eqs. (B.4) and (B.5)]. In particular. 



J5,6(K„r'")=J^r-'^'oosU''-^], 



J.yeixnr'")= J—r-y^cosLry^-^) 
V ^K„ \ 6 / 



Replacing cos(K„r^^^ -27r/3) with sin^ATnr-'^^ -j^jf^^^ can 
make use of the boundary condition 6f(l,t) = to obtain the 
spectrum /(■„ « (n-i-l/6);r, n = 1,2, .... The corresponding 
proper functions are 



Sfnir, t) « exp 



r'l^ 



« + ^j(Mr3/2 + ^(l-At2);j 

- - I ^ . 
cr n y n + 1/6 



Sgnir, t) X exp 



„ + iJ|M,3/2^3^1_;^2)^J 



^-1/2 _ 



n V « + 1/6 



cos 



n + ^ I ^ 



,3/2 , 



(B.8) 



(B.9) 



In Fig. B.l we plot two low radial modes {M - 1/2), that is, 
their real parts, to compare them with their asymptotic expres- 
sions according to Eqs. (B.8)and (B.9). We note that the differ- 
ence between the actual and the asymtotic forms diminishes as 
n or r grow, but it is small even in the most unfavourable case, 
namely, n - I and r < 0.5. 

The simpler form in terms of trigonometric functions al- 
lows us to obtain an analytic form of the orthogonality of eigen- 
functions: 

{yn,ym) ~ i2x 

jgi(m-«);rrJ ^\ ^cos((m - n)nri) - iM sin((m - n)nr^)^dr 



7r2 i(l+6m)(l+6n) (1 - JVP-) 



10 Jose Gaite: Perturbations of self-similar Bondi accretion 

8 

- ^nm /I , ^ X 7 /-I • (B.IO) 

In particular, we have a useful approximation of the norm. We 
have found that the norm given by this formula is in fact a good 
approximation already for n = 1, with a relative error smaller 
than 0.001 (for At = 1/2), which rapidly diminishes as n in- 
creases. 
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